%**********************************************************************
% DISE-ADR: Dynamic Integrated Space Economy model
% DISE-ADR MODEL
% This file: adr.m
% Version 1.0. May 2025
% Paper: "ADR"
% Terminal conditions: Saving rates
% Aneli Bongers and Jose L. Torres (2025)
% Data: Odyssey2010
% Note: CasADi required
% Scenario: Optimal ADR abatement 
%**********************************************************************
clear all
close all
% Importing CasADi space
import casadi.*
% Checking CasADi instalation
x = MX.sym('x')
disp(jacobian(sin(x),x))
clear x;


Params=Odyssey2010;

% List x: Earth capital, Space capital, Derelict satellites, Rocket bodies,
% and Fragments >1cm
% List u: Consumption
% List v: Investment rate on Earth capital, investment rate on Space
% capital, derelict satellites abatement, rocket bodies abatement, and fragments abatement
nx=6;  % Number of states +1
nu=1;  % Number of controls
nv=5;  % Number of endogenous solutions (decision variables)
% Initial values for state variables
x0  =   [Params.k0 Params.s0 Params.W0 Params.Z0 Params.F0 0]'; % Including the sum of discount utility

% Seeds for the decision variables 
vini(1,:)   =   [0.26*ones(1,Params.T+1)];
vini(2,:)   =   [0.001*ones(1,Params.T+1)];
vini(3,:)   =   [0.001*ones(1,Params.T+1)];
vini(4,:)   =   [0.001*ones(1,Params.T+1)];
vini(5,:)   =   [0.05*ones(1,Params.T+1)];

% Compute seed values for trajectories of states and control variables
xini  =   zeros(nx,Params.T+1);
for i=1:Params.T;
    if i==1
        xini(:,1)   =   x0;
    end
    [sa,sb]   =   dynadr(xini(:,i),vini(:,i),i,Params);
    xini(:,i+1) =   sa;
    uini(:,i)   =   sb;
end
uini(:,Params.T+1)    =   sb;
xvuini                =   [xini(:);vini(:);uini(:)];

%*************************************************************************
%Lower and upper bounds
%*************************************************************************
%Decision variables
v_lo=[[zeros(1,Params.T-24) 0.26*ones(1,25)];[zeros(1,Params.T-24) 0.001*ones(1,25)];[zeros(1,Params.T-24) 0.0001*ones(1,25)];[zeros(1,Params.T-24) 0.0001*ones(1,25)];[zeros(1,Params.T-24) 0.01*ones(1,25)]];
%v_lo    =   [zeros(1,Params.T+1);zeros(1,Params.T+1)];
v_up    =   [ones(1,Params.T+1);ones(1,Params.T+1);ones(1,Params.T+1);ones(1,Params.T+1);ones(1,Params.T+1)];
v_lo    =   v_lo(:,1:Params.T+1);
v_up    =   v_up(:,1:Params.T+1);
%State variables
x_lo    =   zeros(nx,Params.T+1);
x_lo(end,:) =   -inf;
x_up    =   inf*ones(nx,Params.T+1);
%Control variables
u_lo    =   zeros(nu,Params.T+1);
u_up    =   inf*ones(nu,Params.T+1);
%All variables
xvu_lo  =   [x_lo(:);v_lo(:);u_lo(:)];
xvu_up  =   [x_up(:);v_up(:);u_up(:)];
%**************************************************************************
% NLP problem
%**************************************************************************
x   =   SX.sym('x',nx,Params.T+1);
u   =   SX.sym('u',nu,Params.T+1);
v   =   SX.sym('v',nv,Params.T+1);
constraints =   SX.sym('constraints',nx+nu,Params.T+1);
for i=1:Params.T
    if i==1
        constraints(1:nx,1) =   x(1:nx,1)-x0(1:nx);
    end
    [sa,sb] =   dynadr(x(1:nx,i),v(:,i),i,Params,u(:,i));
    constraints(1:nx,i+1)   =   x(1:nx,i+1)-sa;
    constraints(nx+1:end,i) =   u(:,i)-sb;
    if i==Params.T
        constraints(nx+1:end,i+1)   =   u(:,i+1);
    end
end
% Social welfare
SWelfare     =   x(end,Params.T+1);

nlp =   struct('x',[x(:);v(:);u(:)],'f',SWelfare,'g',constraints(:));

%IPOPT options
opts                                    =   struct;
opts.ipopt.max_iter                     =   5000;
opts.ipopt.aceptable_tol                =   1e-15;
opts.ipopt.acceptable_obj_change_tol    =   1e-15;
opts.ipopt.mu_strategy                  =   'adaptative';
 
OCP =   nlpsol('OCP','ipopt',nlp);

lbg =   zeros(nx+nu,Params.T+1);
lbg(end-1:end,end)  =   0;
ubg =   zeros(nx+nu,Params.T+1);
ubg(end-1:end,end)  =   0;
sol=OCP('x0',xvuini,'lbx',xvu_lo,'ubx',xvu_up,'lbg',0,'ubg',0);

x_s=full(sol.x);

u_sol=x_s(end-length(u(:))+1:end);
u_sol=reshape(u_sol,nu,Params.T+1);
v_sol=x_s(end-length(u(:))-length(v(:))+1:end-length(u(:)));
v_sol=reshape(v_sol,nv,Params.T+1);
v_sol=v_sol(:,1:end-1);
x_sol=x_s(1:end-length(u(:))-length(v(:)));
x_sol=reshape(x_sol,nx,Params.T+1);

Welfare  =   -x_sol(nx,end);

% Decision variables
srk =   v_sol(1,:);
srs =   v_sol(2,:);
bw  =   v_sol(3,:);
bz  =   v_sol(4,:);
bf  =   v_sol(5,:);

% State variables
k=x_sol(1,:);
s=x_sol(2,:);
W=x_sol(3,:);
Z=x_sol(4,:);
F=x_sol(5,:);

% Control variables
c=u_sol(1,:);

% Rest of variables
for i=1:Params.T
    y(i)=Params.a(i)*k(i)^Params.alpha1*s(i)^Params.alpha2*(Params.N(i))^(1-Params.alpha1-Params.alpha2);
    D(i)=W(i)+Z(i)+F(i);
    S(i)=Params.mu*s(i)/Params.q(i);
    XX(i)=Params.theta*D(i)*S(i);
    L(i)=Params.mu/Params.eta*(1-Params.m(i))*srs(i)*y(i);
    ow(i)=Params.ni1*bw(i)^2*y(i);
    oz(i)=Params.ni2*bz(i)^2*y(i);
    of(i)=Params.ni3*bf(i)^2*y(i);
end

Periods=2023:2222;

% Figures
% Figure 1: Operational satellites
figure;
plot(Periods,S(1:Params.T-50),'LineWidth',2)
grid;
xlabel('Horizon (years)')
ylabel('Satellites')
title('Number of operational satellites')

% Figuere 2: Population of fragments
figure;
plot(Periods,F(1:Params.T-50),'LineWidth',2)
title('Number of fragments >1cm')
figure;
plot(Periods,W(1:Params.T-50))
title('Number of derelict satellites')
figure;
plot(Periods,Z(1:Params.T-50))
title('Number of rocket bodies')
figure;
plot(Params.theta*D(1:Params.T-50))
title('Kessler syndrome')
figure;
plot(Periods,L(1:Params.T-50))
title('Number of launches')
figure;
plot(Periods,XX(1:Params.T-50))
title('Number of destroyed satellites')
figure;
plot(Periods,srk(1:Params.T-50))
title('Investment in Earth capital')
figure;
plot(Periods,srs(1:Params.T-50))
title('Investment in Space capital')

% Save data for Figures
save data5.mat S L D XX W Z F bw bz bf ow oz of y;
